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HIGHLIGHTS 


•  Subbituminous  coal  pyrolysis  was  investigated  using  ReaxFF  molecular  dynamics. 

•  Primary  decomposition  reactions  begin  with  intramolecular  bond  cleavage. 

•  The  formation  mechanisms  for  typical  pyrolysis  products  were  explored. 

•  Phenoxyl  groups  play  a  key  role  in  the  hydrogen  transfer  process  of  gas  generation. 

•  The  reaction  mechanisms  are  coincident  with  previous  experimental  results. 
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A  series  of  molecular  dynamics  simulations  using  the  ReaxFF  reactive  force  field  was  carried  out  to 
investigate  the  mechanism  of  initial  thermal  decomposition  associated  with  pyrolysis  of  a  kind  of 
subbituminous  coal.  The  calculation  results  show  that  the  primary  decomposition  reactions  of  Hatcher 
subbituminous  model  begin  with  intramolecular  changes  such  as  the  cleavage  of  unstable  C— C  and 
C— 0  bonds.  The  formation  mechanisms  for  typical  pyrolysis  products  were  explored.  For  example,  the 
initial  pathway  for  the  formation  of  CO  is  by  the  decarbonylation  of  carbonyl  or  carboxyl  group,  while 
C02  is  mainly  produced  by  hydrogen  transfer  and  decarboxylation  of  carboxyl  groups.  CH4  can  be  formed 
mainly  by  CH3  free  radical  abstracting  a  hydrogen  atom  from  the  hydroxyl  group.  H2  is  formed  by  two 
hydrogen  atoms  from  one  or  two  groups  bonding  together,  which  makes  the  residue  fragments  more  sta¬ 
ble.  Hydrogen  can  also  react  with  oxygen-containing  free  radicals  or  unsaturated  bonds.  Combining  Rea¬ 
xFF  molecular  dynamics  (RMD)  simulation  and  density  functional  theory  (DFT)  calculation,  we  find  that 
the  free  radical  C9H90'  is  an  important  fragment  during  the  pyrolysis  process  of  Hatcher  subbituminous 
model.  As  a  precursor  for  cresol,  it  can  capture  hydrogen  radical  to  form  intermediate  C9Hi0O  and  may 
continue  to  produce  o-cresol  and  ethylene  in  the  presence  of  hydrogen  resource.  These  simulation  results 
for  the  initial  pyrolysis  process  and  the  reaction  mechanisms  agree  with  previous  experimental 
observations. 

©  2014  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Subbituminous  coal,  a  type  of  coal  between  lignite  and  bitumi¬ 
nous  [1],  is  rich  in  volatile  and  has  weak  caking  propensity  so  that 
it  is  highly  suitable  for  pyrolysis  to  co-produce  liquid  (tar),  gas  and 
solid  (char)  products.  Therefore,  there  are  a  lot  of  technical  devel¬ 
opment  works  worldwide  for  establishing  the  processes  of  coal 
pyrolysis  [2-6].  However,  these  technical  activities  are  still  seri¬ 
ously  suffering  from  the  excessive  formation  of  heavy  tar  compo¬ 
nents  which  tend  to  cause  operational  problems  in  pyrolysis  and 
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also  lower  the  liquid  product  quality.  The  solution  to  this  problem 
is  the  control  of  pyrolysis  reactions,  requiring  the  clear  under¬ 
standing  of  the  reaction  process  involved  in  pyrolysis. 

Studies  on  the  mechanism  of  coal  pyrolysis  have  attracted 
much  attention  in  the  past  several  decades  [7,8].  Meanwhile, 
knowledge  of  coal  structure  has  been  made  significant  progress 
and  there  are  more  than  130  coal  models  proposed  since  1940s 
[9  .  The  different  origin  and  different  degrees  of  coalification  result 
in  structural  diversity  of  coal,  especially  for  the  carbon  backbone. 
Several  models  were  published  as  an  aide  to  understanding  the 
pyrolysis  behavior  of  different  rank  coal,  such  as  lignite  by  Tromp 
and  Moulijin  [10  ,  bituminous  coal  by  Solomon  et  al.  [11-13] 
Structural  fragments  were  collected  by  destructive  techniques 
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such  as  flash-pyrolysis  or  wet  chemistry  approaches  that  imposed 
necessary  limitations  on  the  component  molecules.  The  functional 
groups  and  their  changes  during  pyrolysis  have  been  measured 
through  experimental  approaches,  such  as  FTIR  [14,15  ,  NMR 
[16,17]  and  GC-MS  [18,19  ,  and  the  structures  including  heteroat¬ 
oms,  like  sulfur  [20,21],  nitrogen  [22,23]  and  oxygen  [24  ,  have 
been  elucidated  by  X-ray  absorption  near  edge  structure  (XANES) 
spectroscopy.  However,  except  in  situ  IR  for  the  detection  of 
changes  in  the  functional  groups  [14],  there  is  almost  no  adequate 
online  detection  method  to  obtain  more  useful  information  for  the 
continuous  variation  of  species  in  the  process  of  coal  pyrolysis.  In 
addition,  coal  pyrolysis  is  also  a  complex  process  involving  a  large 
number  of  chemical  reactions.  These  chemical  reactions  are  cou¬ 
pled  together,  with  vast  free  radical  intermediates  generated  in 
the  initial  stage.  Such  radicals  are  short-lived  species  so  that  they 
are  hard  to  be  captured  by  traditional  experimental  approaches. 
Instead,  the  computational  method  should  be  capable  of  modeling 
the  chemical  reactions  involved  in  coal  pyrolysis. 

Since  the  reactive  force  field  (ReaxFF)  was  developed  by  van 
Duin  and  Goddard  et  al.  [25-27],  the  computational  capabilities 
for  chemical  reactions  of  large  systems  have  been  significantly 
improved.  The  ReaxFF  reactive  force  field  is  a  general  bond  order 
dependent  potential.  It  has  been  proven  to  be  a  smooth  transition 
from  non-bonded  to  bonded  interactions  and  has  been  demon¬ 
strated  not  only  to  retain  nearly  the  accuracy  of  quantum  chemis¬ 
try  calculations  but  also  to  be  as  low  as  classical  molecular 
dynamics  at  the  computational  cost.  The  method  has  been  success¬ 
fully  used  to  explore  the  combustion  process  of  lignite  and  char 
[28-30],  as  well  as  the  mechanism  of  coal  pyrolysis  [31-35].  For 
the  utility  of  ReaxFF  MD  simulation  on  coal  pyrolysis,  Salmon 
et  al.  [31]  investigated  early  maturation  processes  in  Morwell 
brown  coal  [36]  with  several  functional  models  and  a  macromole¬ 
cule,  respectively.  The  macromolecule  includes  many  oxygenic 
groups,  such  as  methoxyl,  hydroxyl  and  carboxyl  groups.  The  sim¬ 
ulation  reproduced  the  defunctionalization  process  of  aromatic 
side  chains  in  their  offline  experimental  results  [36]  and  showed 
carbon  dioxide  as  the  first  product  formed.  The  effects  of  supercrit¬ 
ical  solvents  (methanol  and  water)  on  coal  pyrolysis  were  investi¬ 
gated  by  Chen  et  al.  [32]  using  a  unimolecular  model  compound 
and  by  Zhang  et  al.  [33]  using  Wiser  [37]  bituminous  coal  model. 
Their  results  showed  that  the  solvent  decomposed  firstly  to  form 
active  radicals,  which  attacked  to  weak  bonds  leading  to  the 
decomposition  of  coal  macromolecular  structure.  More  recently, 
Zheng  et  al.  [34,35]  reported  ReaxFF  molecular  dynamics  simula¬ 
tions  of  bituminous  coal  pyrolysis  using  mixed  models  with  Wiser 
[37],  Shinn  [38]  models  and  other  small  molecules  to  investigate 
the  nascent  decomposition  processes  and  product  profiles.  The 
sequence  of  gas  generation,  the  evolution  of  naphthalene  series 
compounds  and  the  reactions  involving  CH3  and  OH'  were  investi¬ 
gated  to  understand  the  pyrolysis  behavior.  These  works  give  us  a 
useful  insight  for  investigating  the  mechanism  of  coal  pyrolysis, 
but  in  order  to  control  the  pyrolysis  process  and  call  for  high- 
quality  gas  and  tar  products,  more  microscopic  details  of  chemical 
reactions  during  coal  pyrolysis,  especially  for  the  generation  pro¬ 
cess  of  key  gas  (CH4,  H2,  CO,  and  C02)  and  liquid  products  (phenols) 
as  well  as  their  intermediates,  still  need  to  be  further  explored. 
Moreover,  coal  type  is  an  important  factor  which  affects  the  prod¬ 
uct  distribution  of  pyrolysis,  comparing  with  the  structural  feature 
of  lignite  and  bituminous  coal,  subbituminous  coal  contains  less 
ester  and  side  chain  than  that  for  lignite,  and  less  aromatic  rings 
than  that  for  bituminous  coal. 

China  has  a  wealth  of  subbituminous  coal  resources,  accounting 
for  about  one-third  of  total  coal  reserves.  Recently,  we  reported  a 
pyrolysis  technology  to  improve  the  quality  of  tar  using  a  fixed- 
bed  pyrolyzer  enhanced  with  internals  for  a  weak  caking  subbitu¬ 
minous  coal  [39].  In  the  present  work,  to  achieve  the  control  of 


pyrolysis  process  at  the  atomic/molecular  level,  we  implemented 
a  series  of  ReaxFF  molecular  dynamics  (RMD)  simulations  to  inves¬ 
tigate  the  initiation  reaction  pathways  leading  to  the  thermal 
breakdown  of  subbituminous  coal  and  describe  the  formation  pro¬ 
cesses  of  gas  and  liquid  products.  Hatcher  40]  subbituminous  coal 
model  was  selected  for  studying  the  mechanisms  of  pyrolysis.  It 
should  be  noted  that  the  subbituminous  model  was  derived  from 
the  coalification  of  lignin  based  on  13C  NMR  and  pyrolysis  GC- 
MS,  and  the  features  of  this  model  includes  aryl  ether  linkages, 
alkyl  cyclization  at  2,4-position  of  phenol  and  less  side  chain  than 
that  for  lignite  [40].  Although  there  is  no  heteroatoms  N  and  S 
involved  in  this  model,  it  still  well  represents  the  skeleton  for  sub¬ 
bituminous  coal.  In  addition,  the  formation  processes  of  possible 
precursors  for  typical  phenolic  products  were  tracked  from  the 
RMD  trajectories,  and  key  intermediates  as  well  as  transition  states 
were  found  by  assistant  DFT  calculations. 


2.  Simulation  methods 

A  subbituminous  coal  model  reported  by  Hatcher  and  his 
coworkers  [40]  was  used  as  an  initial  model  (Fig.  la).  Firstly,  the 
initial  model  was  built  using  Marvin  program  [41],  and  the  geom¬ 
etry  optimization  of  the  models  was  carried  out  using  Forcite  Mod¬ 
ule  with  the  Dreiding  force  field  in  Materials  Studio  (MS)  package 
[42].  Then  a  minimum  energy  structure  was  used  to  initiate  a  ser¬ 
ies  of  molecular  mechanics  calculations  using  simulated  annealing 
to  generate  several  unimolecular  conformations.  The  procedure 
involved  10  annealing  cycles  from  300  to  1300  I<  under  constant 
number,  volume  and  temperature  conditions  (NVT  ensemble) 
and  geometry  optimizations  were  carried  out  after  each  cycle. 
The  minimum  energy  structure  was  used  in  next  ReaxFF  reactive 
molecular  dynamics  step. 

A  macromolecule  involving  six  molecular  models  with  1110 
atoms  was  immersed  in  a  periodic  box  (30  x  30  x  30  A3) 
(Fig.  lb).  To  achieve  an  appropriate  density  of  subbituminous  coal, 
the  density  of  this  initial  macromolecular  system  was  adjusted  to 
0.78  g/cm3  (26.0  x  25.7  x  26.3  A3)  at  300  K  by  using  NPT  (constant 
number,  pressure  and  temperature)  ensemble  at  pressures  of 
0.2  GPa  with  a  damping  constant  of  500  fs.  Next,  we  performed 
10  ps  no-reaction  MD  simulation  at  300  I<  NVT  ensemble  to  relax 
this  system.  Then,  heat  up  simulations  from  300  K  to  2800  K  at 
rates  of  10,  20, 40  and  80  K/ps  for  the  macromodel  were  performed 
to  determine  the  onset  temperature  of  thermal  decomposition. 
Finally,  ReaxFF  reactive  molecular  dynamics  simulations  were 
performed  using  NVT  ensemble  for  1  ns  with  different  tempera¬ 
tures  at  1000,  1600,  1800,  2000,  2200  and  2400  K,  respectively. 
The  velocity  Verlet  approach  was  adopted  to  integrate  and  update 
Newton’s  equation  of  motion  of  atoms  with  a  time  step  of  0.25  fs.  A 
Berendson  thermostat  with  a  100  fs  damping  constant  was  used 
for  temperature  control.  System  configurations  were  saved  with 
every  0.5  ps  and  separate  molecular  fragments  were  identified  by 
a  ReaxFF  bond  order  upon  0.3.  The  parameters  used  in  the  present 
RMD  study  were  C/H/O  ReaxFF  force  field  parameters  [43  .  The 
energy  conservation  was  evaluated  by  the  NVE  simulations  at 
different  temperatures,  which  is  based  on  the  method  reported 
by  van  Duin  et  al.  [44]  (see  supplementary  material). 

Although  the  time  scale  is  still  many  magnitudes  lower  than 
that  usually  observed  in  experiments,  the  simulation  time  (1  ns) 
is  longer  than  that  performed  in  previous  works  (no  more  than 
250  ps)  for  coal  pyrolysis.  The  temperature  can  affect  the  reaction 
rate  based  on  the  Arrhenius  equation,  especially  for  the  reaction 
with  higher  energy  barrier,  which  will  be  obviously  accelerated 
in  comparison  with  that  with  lower  energy  barrier  [45  .  Thus, 
due  to  the  limitation  on  computational  time  scale,  we  chose  a 
much  higher  temperature  instead  of  experimental  temperature 
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Fig.  1.  Structural  model  for  subbituminous  coal,  (a)  2D  unimolecular  structure  [7  ,  the  circled  letters  indicate  connecting  sites  (e.g.  A— A);  (b)  3D  macromolecule  model  with 
six  molecules  in  a  periodic  box. 


to  simulate  the  pyrolysis  process  of  Hatcher  subbituminous  coal 
model,  which  is  in  the  same  way  as  the  other  pyrolysis  studies 
done  by  using  ReaxFF  reactive  molecular  dynamics. 

All  ReaxFF  reactive  molecular  dynamics  simulations  were 
implemented  in  ADF  software  [46].  The  DFT  calculations  were  per¬ 
formed  by  Gaussian03  program  [47]  using  B3LYP  hybrid  functional 
[48,49]  with  6-3 1G**  basis  set  for  the  geometry  optimization  of  the 
stationary  points.  The  vibrational  frequency  calculations  at  the 
same  level  were  carried  out  to  confirm  each  stationary  point  to 
be  either  a  minimum  or  transition  state  (TS).  Intrinsic  reaction 
coordinate  (IRC)  paths  were  calculated  to  connect  each  TS  to  its 
corresponding  reactant  and  product.  Reported  energies  include  a 
correction  for  zero-point  energies. 

3.  Results  and  discussion 

3.2.  Initial  pyrolysis  process  of  subbituminous  coal 

Heat  up  simulations  from  300  to  2800  K  at  rates  of  10,  20,  40 
and  80  K/ps  were  performed  to  determine  the  onset  temperature 
of  thermal  decomposition  on  the  picosecond  time  scale.  As  shown 
in  Fig.  2,  a  slow  heating  rate  is  beneficial  to  produce  new  fragments 
at  low  temperature.  The  new  fragment  was  firstly  produced  at 
~2400  K  during  the  heat  up  simulation  for  the  rate  of  80  K/ps, 
while  it  was  produced  at  1200  K  with  the  heating  rate  of  10  K/ps. 
From  1200  I<  to  2800  K,  the  slower  the  heating  rate  was,  the  more 


Fig.  2.  Determination  of  initiation  reactive  temperature  using  NVT  RMD  at  heating 
rates  of  10,  20,  40  and  80  K/ps. 


the  number  of  new  fragments  was  when  reaching  the  same 
temperature.  It  was  most  probable  (three-quarters  of  the  chances) 
that  the  initial  pyrolysis  fragment  for  various  heating  rates  was  the 
methyl  free  radical  CH3,  as  we  can  observe  from  the  RMD 
trajectories.  Moreover,  methyl  radical  is  a  critical  precursor  of 
methane,  which  is  an  important  pyrolysis  gas.  For  the  heating  rate 
of  20  K/ps,  the  methyl  free  radical  was  produced  at  1930  K.  It 
should  also  be  noted  that  the  initial  formation  of  methyl  free 
radical  underwent  a  series  of  intra-molecular  bond  breaking 
processes,  as  shown  in  Fig.  3.  This  could  be  ascribed  to  the  intra¬ 
molecular  unstable  tension  of  Hatcher  subbituminous  coal  model 
resulting  from  the  substituent  and  cyclization  at  2,4-position  of 
phenol.  This  structural  unit  is  also  present  in  the  lignite  model  of 
Morwell  coal  by  Nimz  and  Salmon  et  al.  [36,501,  which  is  a  struc¬ 
tural  feature  for  low  rank  coal.  Therefore,  intra-molecular  C— C 
bond  cleavages  occur  in  succession,  and  then  the  free  radical  rear¬ 
rangement  will  happen  among  C563,  C564  and  C631 ,  leading  to  the 
methyl  radical  leaving  and  the  formation  of  a  double  bond  between 
C563  and  C564  (Fig.  3).  The  methyl  radicals,  generated  by  cleav¬ 
ages  of  C— O  bonds  of  methoxy  groups  in  Salmon  et  al.’s  simula¬ 
tions  [31],  were  not  found  due  to  the  absence  of  methoxy  groups 
in  Hatcher  subbituminous  coal  [40  . 

To  investigate  the  kinetics  of  the  thermal  decomposition  pro¬ 
cesses,  we  also  performed  reactive  molecular  dynamics  simula¬ 
tions  using  NVT  ensemble  with  1  ns  and  the  simulation 
temperature  was  at  1200,  1600,  1800,  2000,  2200  and  2400  K, 
respectively,  which  was  chosen  based  on  above  heat  up  simula¬ 
tions.  Fig.  4  gives  a  view  of  the  number  of  fragments  changing  with 
the  simulation  time  under  different  temperatures.  The  number  of 
fragments  is  continuously  rising  when  the  temperature  is  lower 
than  2000  K.  While  the  temperature  is  2400  K,  the  number  of  frag¬ 
ments  rapidly  increases  at  the  beginning  of  RMD  simulation,  then 
decreases  at  about  600  ps,  and  finally  increases  to  98  at  1  ns.  The 
decreases  of  the  number  of  fragments  suggest  that  some  polymer¬ 
ization  reactions  occur  within  the  corresponding  simulation  time. 

After  further  analysis  of  the  product  compounds  observed 
within  the  1  ns  RMD  simulation  at  different  temperatures,  the  evo¬ 
lution  of  product  classes  with  temperature  is  obtained  and  dis¬ 
played  in  Fig.  5.  The  number  of  small  fragments  (Ci— C4  and  non- 
aliphatic  gases)  increased  with  elevating  temperature,  while  the 
number  of  light  tar  fragments  (C5— Ci5  compounds)  increased 
firstly  and  then  kept  being  stable.  The  number  of  heavy  tar  frag¬ 
ments  (C16 — C40  compounds)  decreased  when  the  temperature 
exceeded  2000  K.  The  fragments  of  C5— C40  are  considered  as  tar 
in  this  paper,  as  well  as  shown  in  literatures  [16,34],  which  is  based 
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Fig.  3.  Snapshots  from  ReaxFF  heat  up  MD  simulations  show  the  initiation  of  reactions  leading  to  methyl  free  radical  with  (a  and  b)  and  (c  and  d)  for  intra-molecular  bond 
cleavages  (C562— C566  and  C569— C563),  (e  and  f)  for  generation  of  methyl  free  radical.  The  arrows  located  the  broken  bonds. 


Fig.  4.  Thermal  decomposition  of  subbituminous  coal  model  for  NVT  RMD 
simulations  at  1600,  1800,  2000,  2200  and  2400  K. 


Temperature  (K) 

Fig.  5.  Temperature  evolution  for  products  of  different  classes  obtained  by 
pyrolysis  simulations  using  ReaxFF  RMD  at  1  ns. 


on  their  average  molecular  weight  around  80-500  amu.  Those  spe¬ 
cies  of  boiling  points  below  360  °C  was  considered  as  light  tar  [39], 
corresponding  to  the  carbon  number  of  Q- Cl5  compounds  in  this 
paper.  The  simulation  results  suggest  that  the  pyrolysis  tempera¬ 
ture  is  an  important  factor  to  affect  the  product  composition. 
Low  temperature  favors  tar  species  generation  while  high  temper¬ 
ature  helps  to  produce  small  molecule  fragments.  Therefore,  if  we 
hope  to  produce  tar  by  pyrolysis  of  subbituminous  coal,  the  tem¬ 
perature  of  pyrolysis  should  be  controlled  in  an  appropriate  range, 
not  too  low  or  too  high. 

For  coal  pyrolysis,  it  is  commonly  known  that  a  higher  heating 
rate  has  to  cause  the  higher  yields  of  pyrolysis  oil  and  gas.  This 
appears  inconsistent  with  the  clarification  by  MD  simulations 
shown  in  Fig.  2.  It  in  fact  indicates  that  the  fragments  formed  in 
heating  up  are  not  the  final  experimental  pyrolysis  products  but 
the  primary  pyrolysis  products.  The  former  includes  the  effects  of 
both  free  radical  and  secondary  reactions.  As  one  can  see  from 
Fig.  5,  the  generated  products  in  MD  simulations  are  majorly  gas 
species  below  C4,  while  big-molecule  species  are  not  so  many.  This 
result  shows  that  the  MD  simulations  reproduce  rather  the  charac¬ 
teristics  of  products  without  secondary  reactions  that  are  closely 
related  to  the  field  characteristics  of  flow,  temperature  and  transfer 
in  the  reactor. 

Typical  gas  products  formed  after  1  ns  ReaxFF  RMD  simulations 
of  subbituminous  coal  pyrolysis  under  various  temperature  condi¬ 
tions  were  shown  in  Fig.  6.  Carbon  dioxide  was  formed  at  low 
temperature,  whose  number  increased  firstly  with  temperature 
and  then  fluctuated  in  a  small  range  (~10).  This  is  due  to  that  C02 
is  mainly  derived  from  carboxyl  groups  whose  number  is  a  constant 
for  this  system.  Many  H20  and  CO  were  generated  with  rising  tem¬ 
perature,  especially  when  the  temperature  exceeded  2000  K.  Hydro¬ 
gen  was  observed  in  the  RMD  simulations  at  high  temperatures.  The 
number  of  CH4  also  increased  with  raising  temperature,  and  then 
fluctuated  at  high  temperatures.  It  implies  that  CH4  is  unstable 
and  easy  to  react  with  other  fragments  at  high  temperatures.  The 
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Fig.  6.  Typical  gas  products  generated  after  1  ns  ReaxFF  MD  simulations  for 
subbituminous  coal  pyrolysis  under  various  temperature  conditions. 

trend  in  the  formation  of  these  five  typical  gases  with  the  variation 
of  temperature  is  consistent  with  the  experimental  results  [7,13]. 
These  simulated  results  for  these  gases  at  1800-2400 1<  was 
comparable  with  the  experimental  results  for  a  subbituminous  coal 
at  873-1273  K  [39]. 

3.2.  Major  product  species  analysis 

The  temperature  effect  of  the  chemical  composition  is  shown  in 
Table  1,  which  is  observed  after  the  1  ns  NVT  simulations  of 
Hathcer’s  model  at  1600,  1800,  2000  and  2200  K,  respectively. 
The  decomposition  products  are  divided  into  five  classes  of  com¬ 
pounds,  involving  small  molecular  C4_  compounds,  light  tar  C5_i5 
compounds,  heavy  tar  Ci6_40  compounds,  large  molecular  C40+  char 
and  non-aliphatic  gases.  The  results  show  that  the  number  of  low 
molecular  weight  species  (C15_  compounds)  increased  with  the 
increase  of  simulation  temperatures.  However,  at  high  tempera¬ 
tures  (2000  and  2200  K),  the  large  molecular  weight  species  (Ci6+ 
compounds)  decreased  comparing  with  those  at  low  temperatures, 
indicating  that  large  molecular  species  decomposed  into  small 
molecular  species  at  high  temperatures. 

For  non-aliphatic  and  C4_  compounds,  two  representative  spe¬ 
cies  are  C02  and  CH4.  The  number  of  C02  formation  is  about  10, 
which  is  stable  at  high  temperatures.  There  are  many  C4_  species 
observed  at  the  end  of  the  RMD  simulation  with  the  temperature 
2000  K.  More  interestingly,  C9H90  is  unique  fragment  for  C5_i5 
compounds  after  1  ns  RMD  simulation  with  temperatures  at 
1600  and  1800  K,  suggesting  that  the  C9H90  group  is  liable  to  leave 
from  the  macromolecule.  However,  C9H90  is  not  stable  at  high 
temperatures,  for  example  at  2000  K,  which  can  either  bind  a 
hydrogen  atom  to  form  C9Hi0O  or  react  with  another  group  to  pro¬ 
duce  a  large  fragment.  The  number  of  heavy  tar  composition  of 
Ci 6-40  compounds  increased  from  1600  to  2000  I<  while  decreasing 
at  2200  K,  which  means  that  this  kind  of  compounds  may  either 
decompose  into  smaller  fragments  or  aggregate  with  other  frag¬ 
ments  to  form  larger  fragments.  The  C40+  compounds  decomposed 
into  non-aliphatic  gases  and  tars.  The  weight  for  C40+  compounds 
at  2200  I<  is  greater  than  that  at  2000  K,  indicating  that  these 
compounds  at  2200  K  may  be  formed  by  intermolecular 
polymerization. 

In  order  to  further  study  the  evolution  of  species  during  the 
pyrolysis  process,  the  changes  of  decomposed  compounds  over 
simulation  times  were  analyzed  for  1  ns  NVT  RMD  simulation  at 
2000  K,  as  shown  in  fable  2.  The  number  of  small  fragments 
(non-aliphatic,  C4_  and  C5_16  compounds)  increased  with  prolong¬ 
ing  simulation  time,  while  the  fluctuations  for  the  carbon  content 
of  Ci6_40  compounds  occur  during  last  600  ps.  These  mean  that 


Table  1 

Chemical  molecules  observed  after  1  ns  NVT  simulation  at  1600,  1800,  2000  and 
2200  K  for  Hatcher  model. 


1 600  K 

1800  K 

2000  K 

2200  K 

Non-aliphatic  gases 

4  C02 

9  C02 

11  co2 

10  C02 

1  h2o 

1  h2 

2  CO 

7  CO 

4  H20 

7  H20 

1  h2 

C4_  compounds 

1  ch4 

1  ch3 

1  ch3 

1  ch3 

1  C2H4 

4  CH4 

4  CH4 

8  CH4 

1  CH40 

2  C2H2 

3  C2H2 

1  C2H4 

4  C2H4 

1  c2h3 

1  C2H5 

4  C2H4 

1  CH40 

1  ch2o 

4  C2H20 

1  CH40 

2  C3H6 

2  C2H20 

2  C4H4 

4C3H6 

1  C4H40 

1  c3h2o 

2  C4H40 

C5_j5  compounds 

1  c9h90 

5  C9H90 

1  c6h6 

2  C5H5 

1  c8h6o 

1  c6h5o 

1  C10H10 

1  c9h7 

2  C9H10O 

1  c8h7o 

1  C10H10O 

1  c9h60 

1  c12h10 

1  CioHi002 

1  c12h13o 

1  c„h14o 

1  c15h14o 

1  C13H12 

1  c15h„o2 

1  c14h„o 

1  Q5Hi202 

C j 6-40  compounds 

1  C24H2203 

1  C17H13O3 

1  Ci6Hi203 

1  Ci7Hi302 

i  c25h20o5 

1  C2oH2iO 

1  c19h15o 

1  Ci7Hi403 

1  c26h24o4 

1  c22h17o3 

1  Ci9H1402 

1  Ci9Hi403 

1  c25h23o5 

1  C23Hi903 

1  Ci9H1503 

i  c25h24o3 

1  C3oH3204 

i  c23h22o2 

2  C2oH2203 

1  c34h25 

1  C23H2i  O5 

1  C23H2oO 

1  C32H2503 

1  c34h30o4 

1  C25H1702 

1  C40H33O5 

1  C34H3i04 

1  C25H2i02 

i  c35h27o4 

1  C39H2s04 

C40+  compounds 

1  C6oH54010 

1  C43H37O5 

1  c47h39o6 

1  Cs4Hg40g 

1  C87H82O10 

1  C53H4906 

1  Cg2Hg2Og 

1  CgoHggOg 

1  Ci17Hio60l8 

1  C85H75OH 

1  Ci24Hno015 

1  C9iH830ii 

many  small  fragments  came  from  the  decomposition  of  Ci6_40 
compounds.  Most  of  the  fragments  Ci6-C40  was  unique  due  to 
the  limitation  of  system  scale,  except  for  the  fragments  C20H22O3 
in  Table  2,  whose  number  increased  to  2  at  1  ns,  including  two  iso¬ 
mers  with  the  same  phenol  skeleton  and  different  substituent 
group  at  2-postion  of  phenol,  as  shown  in  supplementary  material. 
Due  to  the  limitation  of  simulation  time,  we  did  not  observed  fur¬ 
ther  variation  of  the  fragment  C20H22O3,  which  may  be  decom¬ 
posed  to  form  phenolic  products  if  increasing  the  simulation 
time.  The  weight  for  C40+  compounds  at  1000  ps  is  greater  than 
that  at  600  ps,  indicating  that  there  are  polymerizations  to  occur 
among  different  compounds.  At  the  beginning  stage  of  RMD  simu¬ 
lation,  more  changes  occurred  between  macromolecules  and  only 
three  kinds  of  small  molecular  fragments  were  produced  at 
lOOps,  including  C02,  CH3  and  C9H90.  The  free  radical  CH3  was 
formed  at  the  beginning  of  MD  simulation,  which  is  a  precursor 
for  producing  CH4.  Although  there  was  no  H2  observed  at  the  end 
of  1  ns  RMD  simulation,  it  would  still  be  present  in  the  middle 
stage  of  this  simulation.  We  focused  on  the  evolution  of  five  com¬ 
pounds  including  C02,  CO  CH4,  H2  and  C9H90,  and  the  time  courses 
of  the  number  of  them  are  shown  in  Fig.  7. 
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Table  2 

Chemical  molecules  observed  with  1  ns  of  NVT  simulation  at  2000  K  for  Hatcher 
model. 


100  ps 

400  ps 

600  ps 

1000  ps 

Non-aliphatic  gases 

2  C02 

7  C02 

9  C02 

11  C02 

3  H20 

2  H20 

4  H20 

1  CO 

1  CO 

2  CO 

1  h2 

1  h2 

C4_  compounds 

1  ch3 

2  CH4 

2  CH4 

1  ch3 

2  C2H4 

1  c2h2 

4  CH4 

1  CH40 

2  C2H4 

2  C2H2 

1  c2h2o 

1  CH40 

4  C2H4 

1  C4H4O 

3  C2H20 

1  C2H5 

1  C4H3 

1  CH40 

1  C4H4O 

4  C2H20 

2  C3H6 

2  C4H4 

1  C4H4O 

C5_75  compounds 

1  c9h9o 

1  C6H8 

1  C6H7 

1  CgHg 

1  C8H60 

1  C8HgO 

1  C8HgO 

3  C9H90 

2  C9H10O 

1  C10H10 

1  C10H10 

1  C10H10 

2  C9H100 

1  C10H10O 

1  C10H10O 

1  C10H10O 

1  c15h13o 

1  C„H10O 

1  Ci2Hio 

1  c15h13o 

1  c12h13o 

1  c15h14o 

1  c15h„02 

CJ6_40  compounds 

1  C18H16O2 

1  C16H1502 

1  CigHi402 

1  CigHi203 

1  C25H21O5 

1  C19Hi402 

1  Ci9Hi402 

1  c19h15o 

1  C25H23O5 

1  C2oHig03 

1  C20H21O3 

1  Ci9H1402 

1  C39H3905 

1  C20H21O3 

1  C20H22O3 

1  Ci9Hi503 

1  C20H22O3 

1  C23Hi803 

2  C20H22O3 

1  C25H21O2 

1  C25H21O2 

1  C23H20O 

1  C24H20O5 

1  C29H2803 

1  C25Hi702 

1  C2gH2304 

1  C35H30O4 

1  C25H21O2 

1  C29H2803 

1  C4oH3805 

1  C35H27O4 

1  C3oH2704 

1  C39H2804 

1  C34H26O3 

C40+  compounds 

1  C4iH3608 

1  C53H47Og 

1  C41H30O5 

1  C47H390g 

1  Cg4H5909 

1  C57H5i06 

1  C41H4i05 

1  Cg2H5208 

1  Cg4Hgi09 

1  C64H530H 

1  C5iH4g05 

1  C7iHg70i2 

1  C87H790i2 

1  C88H790i2 

Fig.  7.  Analyses  of  typical  species  in  ReaxFF  MD  simulations  of  initial  stage  for 
subbituminous  coal  pyrolysis  at  2000  K. 

As  shown  in  the  figure,  the  number  of  C02  increased  rapidly 
before  400  ps  and  tended  to  be  stable  at  about  10  after  600  ps. 
Two  CO  molecules  were  formed  at  140  and  720  ps.  The  number 


of  CH4  increased  with  time  and  stabilized  to  4  at  last.  There  were 
three  H2  formed  during  the  RMD  simulation,  but  they  were  very 
active  in  this  system  and  all  disappeared  at  1000  ps.  The  fragment 
CgHgO  was  formed  before  400  ps  and  then  reacts  with  other  free 
radical  to  produce  large  or  small  compounds.  Two  dominant  iso¬ 
mers  for  CgHgO  appearing  during  RMD  simulations  at  1800  and 
2000  I<  are  shown  in  Fig.  7,  which  will  be  discussed  in  the  following 
text. 

3.3.  Formation  mechanisms  of  typical  gas  and  liquid  products 

The  formation  of  carbon  oxides  (CO  and  C02)  is  mainly  derived 
from  the  dissociation  and  emission  of  oxygenated  functional 
groups  in  the  pyrolysis  process  of  subbituminous  coal.  Previous 
simulation  results  have  proven  that  direct  decarboxylation  con¬ 
tributes  to  the  formation  of  C02  in  both  pyrolysis  and  oxidation 
processes  [29,31  ].  This  formation  process  of  C02  was  also  observed 
in  our  simulations.  For  the  simulation  at  2000  K,  hydrogen  of  the 
carboxyl  groups  firstly  transferred  to  hydrogen  acceptors  that  were 
derived  from  four  types  of  oxygen-containing  groups:  ketone,  phe- 
noxyl  radicals,  enol  radicals  and  alkynol  radicals,  as  shown  in 
Fig.  8a.  Then  11  C02  molecules  would  be  formed  by  decarboxyl¬ 
ation  of  carboxyl  free  radicals.  However,  Behar  and  Hathcer  [51] 
reported  that  the  increase  in  C02  did  not  derive  entirely  from  loss 
of  carboxyl  carbon,  which  also  came  from  loss  of  carbonyl  carbon 
in  NMR  spectra.  Thus,  loss  of  carbonyl  carbon  may  be  firstly  to 
form  the  precursor  of  C02,  for  example  CO,  as  can  be  seen  in  our 
simulations.  Comparing  with  the  scales  of  experimental  time  and 
temperature,  the  RMD  simulation  is  just  suitable  for  studying  the 
initial  process  of  coal  pyrolysis. 

Two  formation  pathways  of  CO  derived  from  different  oxygen¬ 
ated  functional  groups  were  observed  in  our  RMD  simulations, 
which  were  shown  in  Figs.  9  and  10,  respectively.  In  pathway  1, 
CO  derived  from  the  carbonyl  group  was  formed  by  the  cleavage 
of  Cal— CO  bond,  following  the  cleavage  of  Car— CO  bond,  as  shown 
in  Fig.  9.  The  order  of  bond  cleavages  accords  with  the  rule  of  min¬ 
imal  energy  bond  breaking  firstly.  The  bond  energies  for  the  break¬ 
ing  Cal— CO  and  Car— CO  bonds,  using  a  model  compound,  were 
calculated  through  DFT  method  with  the  basis  set  superposition 
error  (BSSE)  correction.  The  bond  energy  for  Car— CO  is  higher  than 
that  for  Cal— CO  bond  (94.93  vs  77.94  kcal/mol). 

The  formation  pathway  2  of  CO  derived  from  the  carboxyl  group 
is  illustrated  in  Fig.  10.  The  carboxyl  group  extracted  a  hydrogen 
atom  from  the  hydrogen  donor  to  form  — C(OH)2  free  radical  group 
(Fig.  10b).  Then  the  ether  linkage  broke  down  to  form  Car— O-  free 
radical,  which  is  a  strong  hydrogen  acceptor  and  obtained  a  hydro¬ 
gen  from  a  hydroxyl  group,  as  shown  in  Fig.  lOc-e.  The  unstable 
— C(OH)2  free  radical  group  lost  a  hydrogen  atom  to  recover  the 
carboxyl  group  (Fig.  lOf).  Next,  a  methyl  free  radical  attacked  the 
carboxyl  group  to  form  methanol  (Fig.  lOg  and  h),  which  is  a  cru¬ 
cial  step  in  the  formation  of  CO  derived  from  the  carboxyl  group. 
Finally,  cleavage  of  Car— C  bond  to  produce  CO.  The  results  for  these 
two  pathways  indicate  that  the  formation  of  CO  may  be  originated 
from  different  oxygenated  functional  groups  and  go  through  a 
series  of  complex  processes,  such  as  C— C  bond  cleavages,  intermo- 
lecular  hydrogen  transfers,  the  ether  bond  cleavage  and  the 
dehydroxylation  of  the  carboxyl  group. 

The  primary  components  of  pyrolysis  gas  contain  methane  and 
hydrogen,  whose  formation  reactions  were  shown  in  Fig.  8b  and  c. 
Three  quarters  of  CH4  molecules  were  produced  by  CH3  free  radical 
abstracting  hydrogen  atom  from  the  phenolic  hydroxyl  group  and 
one  CH4  molecule  was  formed  by  hydrogen  transfer  from  the  alkyl 
group  to  CH3.  Hydrogen  was  generated  by  three  ways:  (1)  in  a 
common  fragment,  two  hydrogen  atoms  of  methyl  directly  leave; 
(2)  two  hydrogen  atoms  are  derived  from  — CH2—  and  — CH3  groups 
in  different  fragments;  (3)  hydroxyl  hydrogen  bound  with  the 
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Fig.  8.  Examples  of  chemical  reactions  for  typical  gases  during  RMD  at  2000  K:  (a-c)  formation  pathways  of  C02,  CH4  and  H2  and  (d)  disappear  reactions  of  H2.  (Rn:  aliphatic 
group;  Arx:  aromatic  group). 


f  :T'< 


Fig.  9.  The  formation  pathway  1  of  CO  derived  from  the  carbonyl  group,  (a)  For  the  intact  carbonyl  group,  (b)  for  the  bond  C241—  C242  breaking  and  (d)  for  the  cleavage  of 
bond  C242— C243  and  generation  of  CO.  The  arrows  located  the  broken  bonds. 


hydrogen  atom  on  vinyl  to  form  H2  and  then  the  electron  rear¬ 
rangement  occurred  to  increase  the  conjugated  double  bonds  in 
the  residual  structure.  As  can  be  seen  from  Fig.  7,  H2  was  formed 
during  the  middle  stage  of  RMD  simulation  at  2000  K,  but  it  disap¬ 
peared  after  1  ns  RMD.  Thus,  we  have  tracked  the  reactions  with 
hydrogen  participation.  Fig.  8d  shows  us  that  H2  reacted  with  oxy¬ 
gen-containing  free  radicals  to  stabilize  them  and  form  hydroxyl  or 
FI20,  which  suggests  that  H2  was  of  high  reactivity  in  this  system. 
During  the  formation  pathway  of  gas  product,  phenolic  hydroxyl 
groups  play  an  important  role  in  hydrogen  transfer,  as  a  good 
hydrogen-donor,  which  can  provide  FT  for  the  formation  of  CH4 
or  H2.  While  a  phenolic  hydroxyl  group  loses  FT  to  become  a  phe- 
noxyl  radical,  it  can  abstract  hydrogen  again  from  the  carboxyl  or 


other  hydrogen-donor  in  the  vicinity,  such  as  the  second  type  of 
C02  formation  reactions  (Fig.  8a).  In  addition,  small  molecules  or 
radicals  have  a  superior  mobility,  which  also  increases  the  chance 
of  collision  with  other  active  fragments,  thereby  improving  their 
reactivity. 

Light  tar  is  the  desired  component  formed  during  coal  pyrolysis. 
The  fragment  C9H90  is  a  pivotal  species  during  the  RMD  simula¬ 
tions  of  pyrolysis  process  for  Hatcher  subbituminous  coal.  We 
observed  that  there  were  several  isomers  for  C9H90‘  free  radical 
during  all  RMD  simulations.  The  formation  pathways  of  two  domi¬ 
nant  isomers  were  studied,  which  may  be  the  precursor  of  cresol 
under  appropriate  conditions.  The  formation  pathway  for  the  iso¬ 
mer  A  of  C9H90  fragment  was  shown  in  Fig.  1 1 .  The  initial  structure 
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Fig.  10.  The  formation  pathway  2  of  CO  derived  from  the  carboxyl  group,  (a-c)  for  the  formation  of  099— H525  bond  and  the  cleavage  of  ether  bond  C26— 033,  (d)  for  the 
transition  state  of  H551  from  a  hydroxyl  group  to  033,  (e  and  f)  for  the  cleavage  of  bond  01 00— HI  84  and  HI  84  binding  with  an  oxygen  radical,  (g  and  h)  for  the  hydroxyl 
leaving  from  the  carboxyl  group  and  the  formation  of  methanol,  and  (i)  for  the  formation  of  CO.  The  arrows  located  the  broken  bonds. 


firstly  went  through  ring  rearrangement  to  make  one  of  the  two  six- 
member  rings  change  to  the  five-member  ring  (Fig.  11  a-c).  Then 
the  C271— C273  bond  broke  down  to  form  a  fragment  containing 
double  rings  with  a  free  electron  at  C273  (Fig.  lid  and  e).  Finally, 
electron  rearrangements  occur  among  the  atoms  of  C273,  C280 
and  C282,  leading  to  the  cleavage  of  C280— C282  bond  to  form  the 
C273=C280  double  bond  and  a  free  radical  at  C282  (Fig.  Ilf). 

Cresol  is  an  important  component  of  tar  produced  from  subbitu- 
minous  coal  pyrolysis  [52,53  .  Thus,  as  a  precursor  of  cresol,  the  for¬ 
mation  pathway  for  the  isomer  B  of  C9H90  fragment  was 
investigated  by  combining  ReaxFF  MD  simulations  with  DFT 
calculations.  As  can  be  seen  from  Fig.  12,  the  radical  C9FI90  with 
2,4-position  cyclization  (RO)  was  observed  to  be  formed  at  the 
initial  stage  of  the  ReaxFF  MD  trajectories,  which  went  through  a 
ring-opening  transition  state  to  produce  the  isomer  (IM1)  of 
C9H90  fragment.  Then  isomerization  will  occur  from  IM1  to  the 
isomer  B  (IM2)  by  a  smooth  transition  state  (TS2)  of  dihedral  twist. 
If  there  are  active  hydrogen  atoms  near  the  reaction  site  in  this  sys¬ 
tem,  IM2  can  abstract  hydrogen  atom  to  saturate  the  2-position 
methylene  radical.  The  active  hydrogen  atoms  may  involve  hydro¬ 
gen  free  radical  or  other  groups  to  easily  lose  a  hydrogen  atom. 
The  isomer  B  of  C9H90  can  bond  with  a  hydrogen  free  radical  existing 
near  the  reaction  site  to  form  a  stable  intermediate  PO,  as  shown  in 


Fig.  12.  The  energy  of  the  transition  state  is  187.5  kcal/mol  based 
on  that  of  C9Hi0O,  calculated  by  DFT  methods.  The  PO  is  one  of  C3 
phenols,  in  accordance  with  the  experimental  data  of  subbitumi- 
nous  coal  pyrolysis  from  Flatcher  [  1 9,40]  and  others  [52,54  .  Because 
the  pyrolysis  process  could  provide  sufficient  energy,  the  4-vinyl 
group  of  intermediate  CgFIioO  (PO)  may  interact  with  5-carbon  in 
aromatic  ring  to  form  intermediate  (IM3)  via  a  transition  state  TS3. 
If  there  are  H2  in  the  system,  one  of  double  bonds  in  IM3  can  be  sat¬ 
urated  via  TS4  to  form  a  low-energy  intermediate  IM4.  Although  the 
energy  for  TS4  is  up  to  114.2  kcal/mol,  it  is  still  lower  than  that  for 
the  reactant  C9FI90  (RO),  which  has  been  observed  in  ReaxFF  MD 
simulations.  Finally  the  intermediate  IM4  can  decompose  into 
o-cresol  and  ethylene  via  TS5.  The  energy  for  products  is 
-4.6  kcal/mol,  lower  than  that  of  CgFIioO,  which  is  a  favorable  pro¬ 
cess  in  energy.  DFT  calculations  only  demonstrated  one  of  probable 
formation  pathways  for  cresol,  but  it  was  a  valuable  approach  for 
understanding  the  formation  mechanisms  of  typical  products. 

This  work  mainly  focused  on  the  formation  process  of  typical 
gas  and  liquid  products  during  pyrolysis.  For  the  NVT  simulations 
at  different  temperatures,  we  thought  that  the  reactions  for  the 
generation  of  pyrolysis  gas  and  liquid  products  can  be  extrapolated 
to  lower  temperatures.  The  scale  of  the  simulation  system  is  much 
less  than  that  for  the  coal  particle,  so  the  reactions  for  large 
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Fig.  11.  The  formation  pathway  for  isomer  A  observed  during  RMD  simulation  at  2000  I<.  Snapshots  (a-c)  show  the  arrangement  of  double  six-member  rings,  including  the 
bond  formations  (C278— C274  and  C279— C275)  and  cleavages  (C278— C279  and  Cl 74— C275).  Snapshots  (d  and  e)  show  the  formation  of  isomer  A  by  the  bond  cleavages 
(C271— C273  and  C280— C282).  The  dashed  lines  indicated  the  bond  formation  and  the  arrows  located  the  bond  breaking. 


Fig.  12.  Energy  diagrams  of  C9H100  formation  and  its  reaction  pathway  in  the  presence  of  H2. 


fragments  (char)  cannot  be  extrapolated  to  lower  temperatures.  At 
high  temperatures,  the  large  fragments  are  liable  to  be  pyrolyzed 
to  form  small  fragments.  Further  simulations  will  be  conducted 
to  compare  the  formation  mechanisms  of  the  same  products  in  dif¬ 
ferent  coal  structures.  In  addition,  the  simulation  system  is  not 
involved  heteroatoms  S  and  N  that  are  often  reactive  in  pyrolysis 
of  most  coals,  which  will  be  the  target  of  future  studies. 

4.  Conclusion 

In  this  article,  we  have  investigated  the  pyrolysis  mechanisms 
of  subbituminous  coal  using  the  ReaxFF  reactive  molecular  dynam¬ 
ics  simulations  assisted  by  DFT  calculations.  The  results  demon¬ 
strate  that  the  thermal  decomposition  is  initiated  by  C— C  bond 
cleavage  following  intramolecular  hydrogen  transfer.  During  1  ns 
reactive  molecular  dynamics  simulations,  there  are  many  frag¬ 
ments  to  be  observed,  including  small  molecule  gases,  precursors 
of  phenols  and  other  species  of  tar.  With  increasing  temperature 


and  simulation  time,  the  second  polymerizations  appeared  among 
different  radical  fragments.  The  initial  pathway  for  the  formation 
of  CO  is  through  the  decarbonylation  of  carbonyl  or  carboxyl 
groups.  Phenoxyl  groups  can  provide  hydrogen  as  hydrogen  donor 
and  lose  FT  to  become  hydrogen  acceptor,  which  play  a  key  role  in 
the  formation  pathway  of  C02,  CH4  and  H2.  The  free  radicals  of 
CgHgO"  are  observed  to  be  an  important  fragment,  as  the  precursor 
of  cresol  and  C3  phenol,  which  can  capture  hydrogen  radical  to 
form  a  relatively  stable  intermediate  C9Hi0O  and  then  continue 
to  produce  o-cresol  and  ethylene  in  the  presence  of  hydrogen 
resource.  The  information  on  mechanisms  and  chemical  events 
during  the  pyrolysis  processes  of  Flather  subbituminous  coal 
model  is  coincident  with  the  results  of  previous  experiments. 
These  results  of  the  study  validated  that  the  approach  of  ReaxFF 
MD  simulation  combined  with  DFT  calculation  can  provide  useful 
insights  into  complicated  reaction  processes  and  give  a  reasonable 
atomistic  description  of  the  initial  pyrolysis  mechanisms  for  subbi¬ 
tuminous  coal. 
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